capture program drop runite
program define runite
	local j=`1'
	local qt=`2'
	use $data\temp012.dta , replace
	if `j'>1 {
		quietly: bsample, strata(time)  // use original sample in 1st loop
	}
	*** calculate RIF ***
	quietly {
		forvalues t=0/2 {
			pctile valx=lrwage3 [aweight=orgwgt] if time==`t', nq(100)  // 100 percentile grids
			kdensity lrwage3 [aweight=orgwgt] if time==`t', at(valx) ///
				gen(evalt denst) width(0.065) nograph	// find density at each grid
			
			if `t'==0 {
				gen rif_`qt'=.
			}
			local qc=`qt'/100
			replace rif_`qt'=evalt[`qt']+`qc'/denst[`qt'] ///
				if lrwage3>=evalt[`qt'] & time==`t'
			replace rif_`qt'=evalt[`qt']-(1-`qc')/denst[`qt'] ///
				if lrwage3<evalt[`qt'] & time==`t'
			
			drop valx evalt denst
		}
	}
	* OB decomposition *** RIF-decomposition ***
	di "doing quantile `qt', rep `j'"
	local qc=`qt'/100
	
	* ssc install oaxaca
	* decomposition without reweighing [E(X_1|t=1)- E(X_0|t=0)]B_0
	local flag=0
	while `flag'!=1 {
		quietly: capture oaxaca rif_`qt' $expvar [aweight=orgwgt] if time==0 | time==1, ///
		detail(groupun:union, ///
		groupoth: nonwhite female  married smsa exper-expe4 edex-ee15 state2-state51, /// 
		grouped: educ2-educ5, groupocc: pubsect manuf partte , ///
		groupind: indu2-indu20) weight(0) swap by(time) relax
		
		if "`e(cmd)'"=="oaxaca" {
			local flag=1
		}
		else{
			quietly: bsample, strata(time)
		}
	}
	
	matrix Ra=e(b)
	di "oaxaca 1/3 completed..."
	
	* composition effects with reweighing [E(X_0|t=1)- E(X_0|t=0)]B_c
	local flag=0
	while `flag'!=1 {
		quietly: capture oaxaca rif_`qt' $expvar [aweight=orgwgt] if time==0 | time==2, ///
		detail(groupun:union, ///
		groupoth: nonwhite female  married smsa exper-expe4 edex-ee15 state2-state51, /// 
		grouped: educ2-educ5, groupocc: pubsect manuf partte , ///
		groupind: indu2-indu20) weight(0) swap by(time) relax
		
		if "`e(cmd)'"=="oaxaca" {
			local flag=1
		}
		else{
			quietly: bsample, strata(time)
		}
	}
	
	matrix Rc=e(b)
	di "oaxaca 2/3 completed..."
	
	* get wage structure effects E(X_1|t=1)*[B_1-B_c]
	local flag=0
	while `flag'!=1 {
		quietly: capture oaxaca rif_`qt' $expvar [aweight=orgwgt] if time==1 | time==2, ///
		detail(groupun:union, ///
		groupoth: nonwhite female  married smsa exper-expe4 edex-ee15 state2-state51, /// 
		grouped: educ2-educ5, groupocc: pubsect manuf partte , ///
		groupind: indu2-indu20) weight(0) swap by(time) relax
		
		if "`e(cmd)'"=="oaxaca" {
			local flag=1
		}
		else{
			quietly: bsample, strata(time)
		}
	}
	matrix Rw=-1*e(b)
	di "oaxaca 3/3 completed..."

	*** collect results ***

	* w/ reweighting
	matrix Rt4=Ra[1,3],Rc[1,4],Rw[1,5],Rc[1,6..10],Rc[1,5],Rw[1,11..16],Rw[1,4], ///
		Rc[1,6]+Rw[1,11],Rc[1,7]+Rw[1,12],Rc[1,8]+Rw[1,13],Rc[1,9]+Rw[1,14], ///
		Rc[1,10]+Rw[1,15],`qc', `j'
	matrix Rt4qtau=(nullmat(Rt4qtau)\Rt4)
	
end